Why?

Randomized quantile residual was proposed in the literature [in 1996] to circumvent the…problems in traditional residuals [for generalized linear models]. However, this approach has not gained deserved awareness and attention, partly due to the lack of extensive empirical studies to investigate its performance.

–Feng et al, 2017

Deep Dive into quantile residuals for model evaluation and the most popular implementation of them, DHARMa.

State of DHARMa

  • DHARMa was first released in 2016. CRAN downloads were sparse until 2019, and began steadily increasing. It averages approximately ~500 downloads per day (source: https://hadley.shinyapps.io/cran-downloads/)

  • It has reverse dependencies with 18 R packaged including glmmTMB, performance, easystats

  • DHARMa does not have a published peer reviewed manuscript paper supporing its release

Standard Residuals from GLMMs Can Be Nonsensical

Linear Mixed Model

\[Y_{ij} = X_i\beta + Z_ja + \epsilon_{ij} \]

\[\epsilon \sim N(0, \sigma^2 \mathbf{I_n}) \] \[a \sim N(0, \sigma_a^2) \]

\[ y_i|a_{j} \sim 𝑁(x_i \beta + a_j, \sigma^2) \]

\[ y_i \sim 𝑁(x_i \beta, \sigma^2) \]

‘Standard’ Residuals

Raw \[\epsilon_i = Y_i - \hat{Y_i} \]

(internally) Studentized:

\[\epsilon_i = \frac {Y_i - \hat{Y_i}} {\hat{s}}\]

Pearson/Scaled:

\[ \epsilon_i = \frac {Y_i - \hat{Y_i}} {sd(Y_i)} \]

GLMM Formulation: Binomial Example

Dependent variable: \(Y_{ij}\) (e.g. survival)

survival: \(p_i = Y_i/N\)

\[Y_{ij}|r_j \sim Binomial(N, \pi_{ij})\]

\[\eta_{ij} = \eta + \alpha_i + r_j\] (\(\eta\) = mean, \(\alpha_i\) = treatment, \(r_j\) = grouping variable)

\[\eta_{ij} = log(\frac{\pi_{ij}} {1 - \pi_{ij}} )\]

Raw/Scaled Residuals Can Lack Diagnostic Capabilities for GLMMs

  • Visual patterns are difficult to interpret
  • Over and under-dispersion are difficult to assess
  • Goodness of fit tests are not valid outside of normal-distriubuted variables


Different distributional assumptions and mathematical mean/variance relationships of some distributions require a different approach

The Original

Screenshot of original paper title

Background

For \(y_1,...,y_n\) responses:

\[y_i \sim \mathcal{P}(\mu_i, \phi)\]

CDF: \[F(y; \mu, \phi)\]

For a continuous \(F\), \(F(y_i; \mu_i, \phi)\) are uniformly distributed, and the quantile residuals are:

\[ r_{i,q} = \Phi^{-1} \left\{F(y_i; \hat{\mu_i}, \hat{\phi} )\right\}\]

Example: Survival Times of Patients

Survival: \(y_i \sim Exp(\mu_i)\)

\[ \text{log }\mu_i = \beta_0 + \beta_1 \text{ log } x_i \]

Quantile residuals:

\[ r_{i} = \Phi^{-1} \left\{ 1-exp(y_i/\hat{\mu_i}) \right\}\]

Discontinuous CDF, F

\[r_i = \Phi^{-1}(u_i)\] \[ u_i \sim Uniform(a_i, b_j]\] \[a_i = lim_{y \rightarrow y_i} F(y; \hat{\mu_i}, \hat{\phi})\]

\[ b_i = F(y_i; \hat{\mu_i}, \hat{\phi}) \]

Method implemented in ‘statmod’

Screenshot of statmod on CRAN

This Method Has Been Around

screenshot of 2003 statmod archive


Randomization is used to produce continuously distributed residuals when the response is discrete or has a discrete component. This means that the quantile residuals will vary from one realization to another for a given data set and fitted model.


–Smythe & Dunn (1996)

DHARMa

(Diagnostics for HierArchical Regression Models) Screenshot of dharma on CRAN

The DHARMa Process

  1. Model something using a generalized linear model with a chosen distribution and link function.

\[\mathbf{Y = X\beta+Za}\] \[ \mathbf{Y|a\sim G(\mu, R)} \] linear predictor: \(\mathbf{\eta = X\beta}\)

fitted value: \(E(Y) = \mu\)

Link function: \(\eta = g(\mu)\)

The DHARMa Process

For each observation in the data set:

  1. Simulate new observations (\(n_{sim} = n_{data}\)) using fitted values as the model distributional parameters (e.g. shape, scale) as appropriate. \[ Y_i \sim G(\hat{\mu_i}, \hat{\phi}) \]

  2. Calculate the quantile for the cumulative distribution function

\[ r_{i} = F(y_i; \hat{\mu_i}, \hat{\phi} )\]

Expectations of the Quantile Residuals

\[ r_{ij} \sim Uniform(0,1) \]

\(r_{ij} = 0\): everything is larger

\(r_{ij} = 1\): everything is smaller

\(r_{ij} = 0.5\): right in the middle

DHARMa runs 250 simulations (per observation) by default, they recommend up to 1000

Poisson Example

\[ (\hat{Y_i} = \lambda = 5,; \quad Y_i = 7 )\]

Quantile Residuals for 1 Observation

Poisson GLMM Example

Create a data set with a random effect (“group”, 10 levels), a covariate (“Environment1”) and a Poisson-distributed response variable (“observedResponse”):

testData = createData(sampleSize = 500)

Analyze correctly and incorrectly:

rightModel <- glmer(observedResponse ~ Environment1 + (1|group) , 
                     family = "poisson", data = testData)

wrongModel <- lmer(observedResponse ~ Environment1 + (1|group) , 
                     data = testData)

‘Standard’ Residuals

Quantile Residuals

Correctly Specified Model

Quantile Residuals

Incorrectly Specified Model

Distributions of Quantile Residuals

Quantile Function of the Standard Normal Distribution

Notes on DHARMa Implementation

  • It depends on the simulation functions build into GLMM package
  • Supported packages: lm(), `glm(), lme4, mgcv, glmmTMB spaMM, GLMMadaptive, brms, & more
  • Commonly, only the last stochastic level (e.g. Poisson) is simulated, conditional on the fitted random effects – basically an omnibus test
  • Residuals can be simulated for individuals predictors and tests can be conducted for individual factors

lme4::glmer() Marginal Model

By default, simulations are conducted on a marginal model:

rightModel <- glmer(observedResponse ~ Environment1 + (1|group) , 
                     family = "poisson", data = testData)
sr <- simulateResiduals(rightModel, re.form = NA, plot = TRUE)

lme4::glmer() Conditional Model

Use re.form argument in glmer() to resimulate the data.

res2 <- simulateResiduals(fittedModel = rightModel, re.form = NULL, plot = TRUE)

glmmTMB Marginal Models

The implementation is a very recent addition and lacks clear documention (see the entire discussion).

data("sleepstudy", package = "lme4")
m1 <- glmmTMB(Reaction ~ Days  + (Days|Subject), sleepstudy)
m1_res <- simulateResiduals(m1, plot = TRUE)

Marginal Model

Set the simulation conditions for a glmmTMB object. The object is altered in place, hence why it’s copied here.

m2 <- m1
set_simcodes(m2$obj, val = "fix")
m2_res <- simulateResiduals(m1, plot = TRUE)

Other Functionality

(Implemented in DHARMa)

  • Uniformity test
  • Quantile location tesy
  • Outlier test
  • Tests for spatial and temporal autocorrelation
  • Test for zero-inflation
  • Over/Under dispersion tests

Best Practices

  • Feng et al (2018) demonstrated for many simulated scenarios that quantile residuals are a better diagnostic tool for checking model distribution, model parameters, over/under dispersion, and overall lack of fit of GLMMs than ‘standard’ residuals.

  • The plot of residuals versus fitted is the most reliable indicator of a properly specified model

  • We do not know their overall sensitivity of this time to modest model misspecifications

  • Results from quantiles tests for uniformity should be treated with caution - no pattern proves a model is correct; likewise, nonrandom matterns are not always a deal-breaker

  • Outliers identified are warnings of possible outliers and lack quantitative information

Sources

  • Bates DM, Maechler M, Bolker BM and S Walker (2015). “Fitting Linear Mixed-Effects Models Using lme4.” Journal of Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.

  • Brooks ME, Kristensen K, van Benthem KJ, Magnusson A, Berg CW, Nielsen A, Skaug HJ, Maechler M and BM Bolker (2017). “glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling.” The R Journal, 9(2), 378-400. doi: 10.32614/RJ-2017-066.

  • Dunn KP, and GK Smyth (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics 5, 1-10.

  • Feng et al (2017). Randomized quantile residuals: an omnibus model diagnostic tool with unified reference distribution. arXiv https://doi.org/10.48550/arXiv.1708.08527

  • Gelman A and J Hill (2007). Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press, New York.

Sources

  • Gerber EAE and BA Craig (2024) Residuals and diagnostics for multinomial regression models. Statistical Analysis and Data Mining: An ASA Data Science Journal 17:e11645. https://doi.org/10.1002/sam.11645
  • Liu D, Zhang H (2018) Residuals and Diagnostics for Ordinal Regression Models: A Surrogate Approach. J Am Stat Assoc 113:845–854. https://doi.org/10.1080/01621459.2017.1292915
  • Pinheiro JC and DM Bates (2000). Mixed-Effects Models in S and S-PLUS. Springer Verlag, New York.

  • Stroup WW, Ptukhina M and J Garai (2024). Generalize Linear Mixed Models: Modern Concepts, Methods and Applications. CRC Press, Boca Raton.

Final Thoughts

  • There is a lack of information on using quantile residuals in a mixed model context and the uses of conditional versus marginal residuals
  • RTFM: read the DHARMa vignette and when needed, function documentation. Not all all R code documentation is well written and helpful, but this package is very helpful.
  • The quantile residuals described herein do not apply to multinomial models. An extension of quantile residuals for these models has been developed, but it lacks an easy implememtation.



(G)LMMs are hard - harder than you may think based on what you may have learned in your second statistics class….


–Ben Bolker, GLMM FAQ

Binomial Example

\[y_i \sim binomial(n, p_i) \]

\[ \text{logit} (p_i) = \beta_0 + \beta_1x_i \] \[r_i = \Phi^{-1}(u_i)\]

\[ u_i \sim Uniform(a_i, b_j] \]

\[a_i = lim_{y \uparrow y_i} F(y; n, \hat{p_i}) =\displaystyle\sum_{i=0}^{\lfloor k-1 \rfloor} \begin{pmatrix} n \\ i \end{pmatrix} p_i(1-p)^{n-i}\]

\[b_i = F(y_i; n, \hat{p_i}) =\displaystyle\sum_{i=0}^{\lfloor k \rfloor} \begin{pmatrix} n \\ i \end{pmatrix} p_i(1-p)^{n-i}\]

Consider Groups

This is the ‘wrong model’ (missing ‘environment’), but the plot looks okay

wrongModel2 <- glmer(observedResponse ~ 1 + (1|group), 
                     family = "poisson", data = testData)
res2 <- simulateResiduals(wrongModel2, plot = TRUE)

Consider Groups

When plotted by the missing group:

plotResiduals(res2, form = testData$Environment1)

Standard Residual Plots